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Abstract 

The present paper is devoted to implementation of the immersed boundary tech- 
nique into the Fourier pseudo-spectral solution of the vorticity-velocity formula- 
tion of the two-dimensional incompressible Navier-Stokes equations. The immersed 
boundary conditions are implemented via direct modification of the convection and 
diffusion terms, and therefore, in contrast to many other similar methods, there 
is not an explicit external forcing function in the present formulation. The desired 
immersed boundary conditions are approximated on some regular grid points, using 
different orders (up to second-order) polynomial extrapolations. At the beginning 
of each timestep, the solenoidal velocities (also satisfying the desired immersed 
boundary conditions), are obtained and fed into a conventional pseudo-spectral 
solver, together with a modified vorticity. The zero-mean pseudo-spectral solution 
is employed, and therefore, the method is applicable to the confined flows with zero 
mean velocity and vorticity, and without mean vorticity dynamics. In comparison to 
the classical Fourier pseudo-spectral solution, the method needs 0(4(1 -|- log N)N) 
more operations for boundary condition settings. Therefore, the computational cost 
of the method, as a whole, is scaled by (iVlog A''). The classical explicit fourth-order 
Runge-Kutta method is used for time integration, and the boundary conditions are 
set at the beginning of each sub-step, in order to increasing the time accuracy. The 
method is applied to some fixed and moving boundary problems, with different or- 
ders of boundary conditions; and in this way, the accuracy and performance of the 
method are investigated and compared with the classical Fourier pseudo-spectral 
solutions. 

Key words: Fourier pseudo-spectral solution; Two-dimensional Navier-Stokes 
equations; Vorticity-velocity formulation; Immersed boundary method; Solenoidal 
velocities; Moving boundary problems 
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1 Introduction 



Fourier pseudo-spectral solution of the vorticity-based formulations of the 
Navier-Stokes equations (NSE) have been used widely in the two-dimensional 
incompressible flow simulations [13]. However, the classical implementations 
are limited to the regular domains with simple coordinate-coinciding bound- 
aries and periodic boundary conditions. Now, recent advances in the immersed 
and embedded boundary techniques have raised hopes of extending the range 
of applicability of these methods to the more general domains and boundary 
conditions [milST] . 

With the best knowledge of the authors, the immersed boundary method was 
applied into the vorticity-stream function formulation of the NSE by Calhoun 
[ll,12j for the first time. In the Calhoun's work the immersed surfaces are 
introduced, and the overall mass balance is satisfied, by definition of an ap- 
propriate distribution of vorticity source term. More or less similar line was 
followed in the work of Russell and Wang [ID]. They decomposed the effects 
of a solid wall into a no-slip condition, satisfied by the Thom's rule, and 
a no-penetration condition, imposed by a boundary element method (which 
satisfied the overall mass balance). In the work of Linnick and Fasel |3l], a 
higher-order compact method was used, and a source term was defined in 
crossing the discontinuities, which was obtained from a jump function. Re- 
cently, Wang et al. [47j applied the direct forcing idea of Mohd-Yusof [37] into 
the vorticity-velocity formulation of the NSE. They added an explicit vorticity 
source term to the vorticity transport equation, which was obtained by tak- 
ing curl of the forcing functions of the momentum equations in the primitive 
variables form of the NSE. However, all the above methods were based upon 
the finite-difference or finite- volume spatial discretization. 
In the pseudo-spectral solutions, the volume penalization is one of the popular 
remedies, which has been used several times for implementation of the no-slip 
condition. The method was first proposed in the primitive variables formu- 
lation of the NSE by Arquis and Caltagirone [4], and then re-formulated by 
Angot [3] . In the next years the method was extended to the vorticity-velocity 
formulation [2911421)28] . and used for the fixed, as well as moving boundary 
problems [3T1I43II29IIT6] . 




A new immersed boundary method is proposed in the present paper, in which 
the arbitrary immersed velocity boundary conditions (including the no-slip 
condition), are introduced into the Fourier pseudo-spectral solution of the 
vorticity-velocity formulation of the NSE, without explicit addition of exter- 
nal forcing functions. Instead of the conventional forcing functions, the im- 
mersed boundaries are implemented by direct modification of the convection 
and diffusion terms of the vorticity transport equation in such a way that can 
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Zero velocity margin (^) 

Fig. 1. The flow domain Oy, together with fixed or moving obstacle(s) Vts, and 
other given- velocity regions f^u, are embedded in the regular solution domain D 
(with regular boundary T/)), via a zero- velocity margin B. 



be implemented to the Fourier pseudo-spectral solutions. 
Particularly, it is aimed to use the zero-mean Fourier pseudo-spectral method. 
Therefore, among many possibilities, our suggested flow conflguration is illus- 
trated in Fig. [Tj The regular solution domain D contains all the fluid-solid 
system, including the moving regions of the fluid with given velocity bound- 
ary conditions (that is f2u), in addition to the flxed or moving immersed solid 
objects. Moreover, note that the fluid-solid system is embedded in D via a 
zero-velocity margin B. This margin is considered to ensure that the flow 
physics can be followed by the zero-mean Fourier series, namely, to ensure 
that the velocity and vorticity flelds remain zero-mean, and the dynamics of 
the mean vorticity is zero, during the time integration. In fact, we found that 
considering such a margin can also be useful in other methods {e.g. the flnite 
difference method), in order to overcoming some difficulties related to the ap- 
propriate vorticity boundary conditions and finding solenoidal velocities (see 
the appendix for more discussions about the role of margin B). In the present 
work, this margin is obtained by windowing of the velocity fields. 
On the other hand, definition of this margin is also in the following of our pre- 
vious work IH] , and the idea of an infinitely flat window function, proposed in 
[8] . However, achieving the spectral rates of convergence is not anticipated for 
the finite Reynolds numbers (which is the subject of this paper); therefore, it 
is not needed infinitely fiat window functions. However, the algebraic rates of 
convergence of the Fourier series are competing in many practical situations. 
Although the shape of the inner boundary of B {i.e. Fg) is arbitrary, in our 
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numerical experiments we used rounded rectangular shapes for simplicity. 
The central core of the method is a conventional pseudo-spectral solver of the 
two-dimensional incompressible NSE in the vorticity-velocity form. At the be- 
ginning if each time step, the desired immersed velocity boundary conditions 
are imposed via direct modification of the velocity fields. Then the modified 
vorticity (which will be called the conditioned vorticity), is constructed di- 
rectly from the modified velocities. The modified vorticity is fed directly into 
the pseudo-spectral solver, which gives the vorticity field in the next time. 
Time integration is performed using the classical explicit fourth-order Runge- 
Kutta method to keep the method as simple as possible, and to avoid some 
difficulties associated with the implicit formulations. In comparison to the 
classical pseudo-spectral method, the present method needs C(4(l + logN)N) 
more operations for boundary condition settings, and therefore, the computa- 
tional cost of the method, as a whole, is scaled by (NlogN). 
The paper is continued by presenting the mathematical formulations of the 
classical Fourier pseudo-spectral method and the suggested algorithm for im- 
posing the immersed velocity boundary conditions. Because of its crucial role, 
the boundary condition setting process is presented in details, in an individual 
section. As our numerical experiments, the method is implemented into some 
fixed as well as moving boundary problems, with the surfaces which are coin- 
ciding and non-coinciding with the regular grids. Finally, some basic questions 
about the validity of the method are addressed in an appendix. 



2 Mathematical Formulation 

The mathematical and numerical frameworks of the method are described in 
this section. Beginning from the classical Fourier pseudo-spectral formulation, 
the suggested modifications for imposing the immersed boundaries, and then 
embedding the solution domain into the regular domain are explained. 



2.1 The Fourier pseudo-spectral formulation 

According to Fig. [l| for a two-dimensional velocity vector u = (ui, U2), defined 
on the regular closure D = {D UTd), the dynamics of the vorticity vector 
u = (0, 0, Uz = uje.z = V X u) is obtainable from time integration of the 
vorticity transport equation 




for x G -D, 



in D X (0, T] 



(1) 
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while the velocity vector u satisfies the following Poisson's problem with 
Dirichlet boundary conditions: 



V^u = X Vw in 



(2) 



One of the advantages of the above vorticity-velocity formulation, in compar- 
ison to e.g. many primitive variable formulations, is the possibility of decom- 
position of the kinematics and dynamics of the fiow field at each time instant. 
In fact, for any arbitrary distribution u G C^{D), the physical (divergence- 
free) velocity vector is obtainable from solution of Eq. (|2|, if the appropriate 
boundary conditions are imposed (see the appendix or [20] for more discus- 
sions). As it will be seen, this issue has a vital role in construction of the 
physical immersed velocities in the present method. 

On the other hand, to improve the efficiency of the computations, it is conve- 
nient to change the vorticity transport equation ([T]) to 

d,u = uV^- Q^y. - + - ^)"^"2' (3) 

L ^ ' V ^ / 

Ni Na 



which in comparison to the classical form ([T|, saves one fast Fourier trans- 
form (FFT) in the pseudo-spectral algorithm. Although this formulation has 
been used in some other studies, to the best knowledge of the authors, in the 
pseudo-spectral solution of the incompressible fiow, it was proposed by Chas- 
nov [I5] for the first time. The diffusion term L and the non-linear terms, Ni 
and N2, are named for the future references. In fact, the immersed velocity 
boundary conditions will be introduced to the solution by direct modification 
of these terms. 

For the periodic boundary conditions, the Fourier series provides such an ac- 
curate and efficient tool which makes it worthwhile re-formulating the problem 
in the Fourier space. In this way, the vorticity transport equation dsj) recasts 



dtCj = — z/|kp(u) — kik2{ui — U2) + {k^ — kl)uiU2, 

t Ni N2 

uj{k, t = 0) = uq, 



while the Poisson's problem ^ simplifies to 
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In these equations, (■) stands for the quantities in the Fourier space, k = 
{ki,k2) is the wavenumber vector with the magnitude |kp = kf + k\\ while 
k"*- = {—k2,ki) is the perpendicular wavenumber vector, and i"^ = —1. Practi- 
cally, in the finite dimensional calculations, the nonlinear terms Ni and N2, 
are constructed (and are de-aliased) in the physical space — the algorithm 
which is known (and will be referred in this paper) as the pseudo-spectral 
method. Discretization in time, and time integration of the fully discretized 
system can be done using an appropriate time marching method. The explicit 
fourth-order Runge-Kutta method is used in the present work. 
More or less similar formulations have been used in many efficient and ac- 
curate solvers of the periodic fiow in regular domains. In the sequel we will 
modify equations ^ and ([s]), and suggest an algorithm to use them in the 
fiow configuration of Fig. [T] 



2.2 Implementation of the immersed velocity boundary conditions 



In the present method, without explicit addition of a forcing function in the 
right hand side of Eq. Q, the immersed surfaces are introduced by modi- 
fication of the Ni, N2, and L. Particularly, it is desired to carry out these 
modifications such that the velocities remain solenoidal, and in a manner that 
the method can be implemented easily into a Fourier pseudo-spectral solver. 
The suggested procedure is summarized in Fig. |2} in which the boundary con- 
ditions box (BC) can be more explained by the following remarks: 

1. Given the vorticity field a)", from the initial condition or the last timestep, 
the velocity vector u in the regular domain D, is obtained from Eq. ([s]) and 
two inverse FFTs (box (X)). 

2. The velocity vector u is modified to satisfy all needed immersed velocity 
boundary conditions (box (XX)). This conditioned velocity will be called 

The modifications are carried out by local extrapolations, extension, and 



windowing of u; and will be explained in details in § 2.3 Note that u^*-^ is 
neither necessarily solenoidal, nor its mean value is necessarily zero at this 
point. 

3. The conditioned vorticity w^*^ is re-calculated from u^*-^ (that is, w^*^ = 
V X u^*-'), as it is shown in box (XXX). 

There are two main reasons for this step. Firstly, the solenoidal velocities can 
be obtained from this conditioned vorticity in the next steps, provided that 
the appropriate boundary conditions are implemented (see the appendix 
and the discussions in there); and secondly, the vorticity will be needed in 
the subsequent pseudo-spectral algorithm. 

Although u^'~" can be found by any method {e.g. the finite difference), to 
preserve the spectral accuracy, and because the vorticity in the Fourier 
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Eq. (5) 

(ui,U2] 



2 X IFFT 



(Mi,M2j 



Immersed boundaries 
setting 



{ITT) 



U 

1 



BC 



2 X FFT 



Eq. (6) 



/3C 



Pseudo-spectral solver 



(PS) 



Boundary conditions 

(HC) 



Fig. 2. The main steps of the proposed algorithm. The box (PS) contains a clas- 
sical Fourier pseudo-spectral solver (calculation of the right hand side of Eq. (Q, 
solution of Eq. ([s]), de-aliasing, time integration. . . ). Therefore, a classical Fourier 
pseudo-spectral solution with periodic boundary conditions can be retrieved by 
switching-off the boundary condition setting box, and following the dashed line. 



space is needed in the subsequent steps of the pseudo-spectral algorithm, 
calculation in the Fourier space is suggested here. In this way 

Note that it is aimed to simulate the confined flows, and therefore, the vor- 
ticity fleld has zero-mean according to the Stokes theorem; the fact that le- 
gitimates use of the above equation. Moreover, note that w^*" = IFFT{c2;B*-'} 
is automatically deflned on I), it is double periodic, and it has zero mean — 
the properties that makes it ready to use for the subsequent Fourier pseudo- 
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spectral steps. 

4. The conditioned vorticity u^'" is fed into the classical pseudo-spectral pro- 
cedure (that is, box (PS) in Fig. |2]), in which the solenoidal velocity vector 
ufoi is calculated from 

and these velocities, in addition to the conditioned vorticity are sub- 
stituted into the modified vorticity transport equation 

— BC ~_BC 

c:,(k,t = 0) = c2;^c^ 

Time integration of the above equation yields the new vorticity field w""*"^ 
which closes the algorithm loop. 

An attractive feature of the above procedure is that the boundary conditions 
box can be added easily to a classical pseudo-spectral solver without any 
change in the internal structure. In fact, the box (PS) contains all steps of 



a classical pseudo-spectral solution (as it was formulated in ^2.1). Moreover, 
note that in addition to the solid obstacles, the given immersed fluid velocities 
can be implemented as well. 

On the other hand, as it was mentioned earlier, implementation of the bound- 
ary conditions (that is, box (XX)), can imposes some discontinuities to the 
velocity field. Although the next steps will remove these discontinuities (as it 
will be proven in the appendix), the boundary conditions will change a bit. 
Similar problem was observed in many other immersed boundary methods, 
and caused emergence of some methods like the multi-direct forcing method 
|Tril35] . or the implicit forcing method [3211^ in the last years. In the present 
work, we repeat the above algorithm Mr times to overcome this difficulty; 
where 1 < Mr < 3 depending on the flow. By Mr times repetition of the 
boundary condition settings (and also by development of the solution), the 
velocity field converges to a solenoidal field which satisfies (approximately) 
the desired velocity boundary conditions. 

In our numerical experiments that are presented in this paper, the fourth-order 
Runge-Kutta method is used for time integration, and in order to increasing 
the time accuracy, the boundary conditions are set at the beginning of each 
Runge-Kutta sub-step. 



2.3 Boundary conditions setting 



This section is devoted to a full description of the method used in setting 
of the immersed velocity boundary conditions (that is, modifying u into u^*-^. 
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Fig. 3. The numerical boundary Fn corresponding to the physical boundary Fp. 
The immersed boundary conditions Upb are given on Fp. 



mentioned in item 2 of section 2.2 ). The method includes a local extrapolation, 
borrowed from the finite difference-based immersed boundary methods [H], 
followed by an extension and windowing process, which we used in our pre- 
vious work IH] . Particularly, this combination is chosen in order to achieving 
different (desired) orders of local accuracies in boundary condition settings, 
and the needed flexibility in treating the moving boundary problems. In this 
way, the process of boundary condition setting is divided into the following 
sub-steps: 

(1) Identification of numerical boundary points. 

(2) Evaluation of the numerical boundary conditions. 

(3) Embedding in the regular domain D. 

The details of the above sub-steps are in order. In what follows, we will consider 
one immersed body. For the multi-object problems, the method can be applied 
exactly in the same way. Moreover, according to Fig. |3} we assume that the 
physical velocity boundary conditions Upb are given on the physical boundary 
Fp, and the solution is sought in D \ {QU B), and the regular domain D is 
overlaid by a uniform Cartesian grid {xi,yj). 



2. 3. 1 Identification of numerical boundary points 

In the present method, all calculations are performed on a fixed Cartesian grid. 
Therefore, in the following of our previous work [H], we define the numerical 
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Fig. 4. A numerical boundary point defined on a uniform grid. 

boundary points, which play the role of Eulerian points in some fluid-solid 
interaction methods, which use both the Eulerian and Lagrangian points (see 
t-g- [46]). 

Definition 1. A numerical boundary point corresponding to the given phys- 
ical boundary Fp, is a point {xi,yj) in the Cartesian grid, if and only if 

i) ixi,yj) e{n = nuTp) 

ii) Cij contains at least one point from D\Q 

where Cij is a circle of radius = min(Ax, Ay), centered in {xi,yj). The 
definition is illustrated in Fig. |4] (for a uniform grid), and for some more de- 
tails one can see [S]. The set of all numerical boundary points will be called 
the numerical boundary Fn, and that part of the Cartesian grid which is sur- 
rounded by Fn will be called the numerical immersed domain ^n. 
For the fixed boundary problems, it is just needed to determine the numerical 
boundaries once for all computations, while for the moving boundary problems 
they should be updated in the beginning of each timestep with the computa- 
tional cost of 0{Mq), where Mq is the number of grid points in a box, contains 
the immersed domain = f^Ar U Fjy. 

2.3.2 Approximating the boundary conditions 

Given the numerical boundary points Fn, the next step is to evaluating ve- 
locities on these points (will be called u^^), such that the physical boundary 
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Fig. 5. The geometric parameters, and the interpolation points for a second-order 
boundary condition setting (i.e., Afp = 2). The auxiliary quantities and u^^ are 
obtained from linear interpolations between the neighboring points, which are shown 
by X. 



conditions Upb be satisfied approximately. Among several possibilities, an ex- 
trapolation method is borrowed from the finite difference-based immersed 
boundary methods because of its flexibility, and its consistency with our 
next steps. 

Basically, the method is an order A/^ polynomial extrapolation, along the local 
normal direction n. According to Fig. [sj given the velocity vector u = U2), 
the auxiliary velocities u\ u^^. . . (number of them depends on A/^), are deter- 
mined from some linear interpolations between the neighboring points on the 
Cartesian grid. Then the desired approximation is obtained from 

Ufn =fA/'p(u\u",...;Upb), (9) 

where is the polynomial extrapolation function of order Mp , and Upb is the 
given physical velocity boundary condition. The neighboring points which are 
involved in calculation of the auxiliary velocities are depended on the slope 
of the normal direction n (see In the present paper A/^ = 0,1,2 are 

used in our different test cases, and also in different regions of each test case, 
depending on the desired local accuracies. 

As an example, for a second order extrapolation (that is, A/'p = 2), according 
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to Fig. [5} the numerical velocity boundary conditions can be obtained from 



(5l + 52)P 



Upb + 



53(^1 + S2) 



U — 



(10) 



where S = 61 + 62 + S^', and u , and u are found from linear interpolations. 
With regard to these extrapolations, the following practical points are worth 
mentioning: 

(1) Performing these extrapolations showed to be not so easy for complex ge- 
ometries and moving boundaries. Moreover, for Afp > it may increase the 
sharpness of the velocity fields, which particularly for the spectral solutions 
may increase the Gibbs oscillations by triggering the higher wavenumber 
modes. 

(2) For Mp = 0, it is not any extrapolation indeed, and the boundary condition 
setting reduces to using a mask function, which substantially simplifies the 
process; and in many cases, reduces the Gibbs oscillations. However, for 
the boundaries which are not coinciding with the Cartesian grid, it means 
zero-order implementation of the boundary conditions, which reduces the 
global accuracy of the solution, as it will be seen in § [3j 

By substitution of upj^ into u, an auxiliary modified velocity vector Unb ob- 
tains, which is identical to u in D \ Q-^, and satisfies (approximately) the 
physical boundary conditions on Tp. 



2.3.3 Embedding in the regular domain 

To take advantages of the pseudo-spectral solutions, the auxiliary modified 
velocity vector Unb should be embedded in the regular domain Z); and in order 
to minimizing the Gibbs oscillations, the final modified embedded velocities 
u^^ should be as smooth as possible. Therefore, referring to Fig. [sj the final 
modified velocity u^^ should satisfy the following conditions 



Moreover, smooth transitions from Unb to Upb (in I^n), and from Unb to zero (in 
B) are desired, in order to minimizing the Gibbs oscillations. In this context, 
the following two-step procedure is followed for the sake of flexibility: 

(1) Extension: Given the extrapolation functions fyv^ (for different points 
of Fn and Fg), the desired extensions of Unb (in VLn and are eas- 





(11) 



in B. 
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Fig. 6. A sectional view of an extended and periodized field for J\fp = 1. The are 
obtained from multiplication of a linear extension (with a slope which is computed 
from the internal points adjacent to the physical boundary Tp) by a shifted error 
function. 

ily obtainable. Note that these extensions have C-^^ smoothness (in the 
directions normal to the immersed boundaries), because they are poly- 
nomials of order Mp. The extended velocity will be called (since it is 
defined on D), and note that is not necessarily periodic. 
(2) Smoothing and periodization: Finally u^*-^ is obtained from smooth- 
ing and periodization of u"^. Among some methods that can be seen 
in the literature (see e.g. j^IOl'Q]), we used multiplication of by an 
appropriate window function [41j 



where the window function ^(x) is constructed such that it is sufficiently 
smooth in VL^ and is zero in B. To construct the desired window func- 
tion, we used shifted and scaled error functions, exactly the same as our 
previous work [41j. 

A sectional view of an extended and periodized velocity field for linear extrap- 
olation {Mp = 1) is shown in Fig. |6| The physical boundary condition Upb was 
given on the physical boundary Tp ; and note that because of the slope of the 
solution in the vicinity of Fp (since Afp > 0), a margin with u ^ Upb is ap- 
peared in the immersed body, which will be ignored in the solution procedure, 
like other immersed boundary methods. 

Apparently, the computational cost of boundary condition setting is of 0{Affi'), 




(12) 



13 



where Mqi is the number of grid points in Q' = Q-^ U B. However, It should 
be noted that for the boundaries which are coinciding with the Cartesian grid 
{e.g. the rectangular boundaries for the Cartesian grids), the extrapolation 
step can be bypassed. 



2.4 Filtering of the oscillations 



Although the modified velocities u^^ are C-'^p in the directions normal to 
the immersed boundaries, and they are periodized; however, for the finite 
Reynolds numbers and complex geometries there is no any guarantees for their 
global smoothness [H]. Moreover, the conditioned vorticity cu^*" = V x u^^ 
is generally sharper than u^^. Therefore, as our numerical experiments have 
confirmed as well, the solution suffers from the Gibbs oscillations almost in 
all practical situations. Fortunately, there are some evidence which show that 
despite these oscillations the solutions are still accurate in the weak sense, 
and higher orders of pointwise accuracies can be recovered (if desired), using 
appropriate numerical filters P7|28f31j . 

In the present paper, the Helmholtz filter is used (among many kinds of the 
numerical filters), because of its simplicity, and since it is aimed to use it 
just as a postprocessor in presentation of the results, not for achieving the 
pointwise accuracies or improvement of the rates of convergence. However, 
our other numerical experiments (are not presented in this paper, because 
of lack of any theoretical support), have shown that use of this filter in the 
solution procedure may reduces the Gibbs oscillations. 

Given the non-filtered quantity 0(x), the Helmholtz-filtered quantity 0(x) is 
obtained from the following convolution 

0(x) = (l-«2^2)-i^(^)^ (^3) 



where a is the free parameter of the filter, standing for an appropriate length 
scale (see |25], or the series of works on the a modeling of turbulent flow e.g. 
121]). Since the Laplacian operator is diagonal in the Fourier space, the above 
convolution takes a simple form in the wavenumber space 



where k is the wavenumber vector. Since the dimension of a is length, it is 
convenient to define a more general non-dimensional filtering factor Ca such 
that 
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Fig. 7. Geometric parameters of the dipole-solid wall collision problem. In order to 
maximizing the active grid points, 5 is chosen as small as possible, and different 
smoothness of windowing is obtained by choosing different A' (see also Fig. [s]). 

where L is the length of the solution domain, and 



is the maximum wavenumber, involved in the solution. 

The above filter is used in the next section, wherever it is needed to remove 
the Gibbs oscillations, for better presentation of the results. 



3 Numerical experiments 

In order to assess the capabilities of the method, three test cases are analyzed 
in this section, including the fixed and moving boundaries, as well as given 
non-zero immersed velocity boundary conditions. 

3.1 Dipole-solid wall collision 

As a coordinate-coinciding immersed boundary problem, the recently re-analyzed 
dipole-solid wall collision problem is chosen as our first numerical ex- 

periment. In a series of papers, the problem has been analyzed numerically for 
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Fig. 8. A sectional view of two window functions: a fairly sharp window ^eC^) that 
A = = 6, and a moderately smooth one Qi2{^) that A = 12. The rising distance 
A' is the distance between the positions that g = and g = 1 with the tolerance of 
0(10-14). 



the normal and oblique collisions and a fairly wide range of Reynolds numbers; 
via finite difference, Fourier and Cliebyshev spectral methods, and the volume 
penalization of the NSE [TU|18lll9|ll7f28j . In order to make quantitative com- 
parisons Re = = 1000 is chosen in this paper, similar to the Keetles et 
al. test case l28l . 



3.1.1 The problem setup 

The geometric parameters are illustrated in Fig. [7} where it is chosen 2W = 2 
for simplicity (c.f. |28]). The margin 6 should be chosen as small as possible, 
in order to maximizing the number of active grid points. On the other hand, 
it should be sufficiently wide in order to satisfy the requirements that are 
discussed in the appendix. Here 6 = lOAx is chosen (after a trial and error 
process) in all runs. Although fairly acceptable solutions were obtained for 
smaller 6, solutions exhibited more Gibbs oscillations. 

As it was mentioned in §2.3.3 for the coordinate-coinciding boundaries, the 



extrapolation step can be bypassed, and therefore, the boundary condition 
setting is reduced to merely extension and windowing. In the present test 
case N'p = 1 is chosen; and various numerical experiments were performed 
on a fairly wide range of smoothness of the window functions. The results 
of two typical ones, that is, a sharp window with A = ^ = 6 (henceforth 
will be called QQix.)), and a moderately smooth one A = 12 (which will be 
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Fig. 9. Scatter diagram of the window function gi2{x). The exponential rate of 
decaying, especially in the higher wavenumber, is noticeable (c.f. [8]). 

called (x)), are reported here. A sectional view of these window functions 
are shown in Fig. |8| It is worth mentioning that although a square, with 
sharp corners, is not an ideal shape for a window function; sufficiently smooth 
square-shaped window functions are obtainable. To emphasize this issue, the 
scatter diagram of the absolute values of the Fourier coefficients for ^12 (x) is 
shown in Fig. |9j As one can see, the exponential rate of decaying is observable 
in the higher wavenumber, which can be interpreted as evidence of smoothness 
of the window function |8]. 

In the following of [2H] , the initial condition was constructed by summation of 
two identical mono-poles 

a;o = We[l-(-)']exp[-(-)2], (17) 

placed at {{xi,yi), (x2, 1/2)} = {(0, 0.1), (0, —0.1)}, in which tq = 0.1, and r = 
y/x"^ + y"^. Furthermore, to make quantitative comparisons, the total kinetic 
energy E(t), and the total enstrophy Z(t) are defined as 

E(t) = ^J u2(x, t)dn; Z{t) = ^J a;^(x, t)dn. (18) 
n n 

The main physical parameters of the initial condition are given in Tab. [TJ 
The solutions were performed on a 512^-points uniform grid and the calcu- 
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Table 1 

The physical parameters of the initial condition of the dipole-solid wall collision 
problem. 





V 


u 


E 


Z 


Re 


299.528385375226 


0.001 


1.0 


2.0 


800 


1000 



lations were de-aliased using the |A^-rule. The explicit fourth-order Runge- 
Kutta method with a constant timestep At = 10~^ was used for time integra- 
tion. By definition of the CFL number as 

CFL=||u|U-^, (19) 



in which ||u||oo is the maximum velocity magnitude at each time instant, our 
calculations showed that CFL < 0.5 for all times, and all runs. 



3.1.2 The captured physics 

The contour plots of vorticity fields at a number of time instances are presented 



in Fig. 10 for ^^{'x.) windowing. The time instances are chosen similar to Kee- 
tels et al. |28] for comparisons. Obviously, the no-slip and the no-penetration 
conditions, and all phenomena of the dipole-solid wall collision are observable 
(c.f. [28] or [T7]). More quantitative comparisons can be made by observa- 



tion of time histories of the total kinetic energy and enstrophy. In Fig. 11 
these quantities are presented together with the results of Fourier-Chebyshev 
method and the volume penalization method (these results are extracted di- 
rectly from [28J). As one can see in the left panel, decay of the total kinetic 
energy shows a small deviation from the Chebyshev-Fourier results after the 
first collision. However, this deviation is compensated in the proceeding times, 
such that after about t = 1, it is vanished. In the other words, the dissipation 
rates have been over-predicted during the first collision, and under-predicted 
for the proceeding times (in comparison to the Fourier-Chebyshev solution). In 
our opinion, these differences can be interpreted with regard to the smoothness 
of the solution (which is mainly depended on the grid resolution and sharpness 
of the windowing, for a fixed Reynolds number), and the order of implementa- 
tion of the boundary conditions Afp. In fact, the first extra-dissipation in the 
collision time interval is supposed to be related to appearance of the Gibbs os- 
cillations, which transfers a part of kinetic energy into the higher wavenumber. 
It is this high-wavenumber energy part which will be dissipated, faster than 
it should, in the proceeding times. On the other hand, the under-predicted 
dissipation rates in the next times is presumed to be a consequence of smooth 
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(d):t=0.6 (e);t=0.8 (f):t=1.2 



Fig. 10. Contour plots of the vorticity field for Re=1000, obtained from ap- 
plication of the proposed method on a uniform 512^-point grid with constant 
timestep At = 10~^ (amounts to CFL<0.5 for all times). The contour levels 
(-270,. ..,-50, -30, -10, 10, 30, 50,. ..,270), and the time instants are chosen similar to Kee- 
tels et al. [281. 



walls, instead of sharp boundaries. Unfortunately, direct verification of these 
speculations found to be so difficult in practice, because the smoothness of the 
solution and the order of boundary conditions are not completely independent 



see discussions in § 2.3.2 for some more details). 



In addition to the results of Fourier-Chebyshev method, the results of penal- 
ized NSE are presented in Fig. 11 as well. Obviously the penalization method 
showed more dissipations than our method. Again we suppose that it is a con- 
sequence of sharpness of the generated boundary layers in the penalized NSE 
(with 0{^/eu) thickness [Tif28] ). which is sharper than the window functions 
in the present method (with maximum sharpness of 0{N~^)). Therefore, it 
is not surprising that the penalization method is more dissipative than the 
present method in general. 

On the other hand, in calculation of the total enstrophy (see the right panel 
of Fig. 11 ), this sharpness of the penalized NSE, helped the method in gener- 
ation of the maximum possible vorticity in the Fourier spectral method (that 
is, the mean value of the vorticity in the vicinity of the discontinuities, as 
it is argued in [2B]); while smoothness of the present method resulted in an 
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time(Sec) time (Sec) 

Fig. 11. Comparison of the total kinetic energy E(t) and enstrophy Z(t) of the 
present method with the Chebyshev-Fourier method and the volume penalized NSE 




under-predicted vorticity generation. 

For more verification of the effects of smoothness of the windowing process, 
time histories of the total energy and enstrophy for Qelx) and ^12 (x) are com- 
pared with the penahzation method in Fig. 12 Obviously the dissipation rate 
has changed slightly, while the vorticity generation has affected drastically, 
especially in the collision times. 

As a conclusion, our results showed that for this grid-coinciding immersed 
boundaries, sharper window functions yielded better implementation of the 
boundary conditions, but more dissipation rates. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

time(Sec) time (Sec) 



Fig. 12. Time histories of the kinetic energy and total enstrophy for two window 
functions (with different sharpness Qe{x.) and £)i2(x)), are compared with the results 
of the penalization method. The sharper window function qg{x.) resulted in more 
dissipation and bigger total enstrophy than the smooth window function ^>i2(x). 
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Fig. 13. The rates of grid convergence of the vorticity field for three different window 
functions at t = 0.35. The rate of convergence approaches to N~'^ as the sharpness 
of the window function increase. 



3.1.3 Spatial and temporal rates of convergence 

It is a well-known fact that the rate of convergence of a spectral solution is 
highly depended on the smoothness of the solution. Particularly, for the non- 
smooth fields that the Fourier series are suffered from the Gibbs oscillations, 
the £^ rate of convergence reduces to about 0{1). Fortunately, despite this 
poor rate of convergence, it has been shown that high orders of pointwise 
convergence can be recovered by use of some appropriate numerical filters 
[271l28f31j : although the rates of convergence are depended on the order of 
the numerical filter this time (confirming that loss of uniform convergence 
does not mean necessarily loss of the accuracy in the spectral solutions). In 
this section, to observe the rate of convergence of the method separately, the 
results are reported without any kind of filtering. However, like some other 
spectral methods j28l|3T] . these rates can be improved using the numerical 
filters. 

For the initial times that the effects of the solid walls are negligible, the 
exponential rates for the grid convergence are easily observable for the regions 
far from the walls. Therefore, merely the rates of grid convergence in the 
collision time [t = 0.35 in particular) in the near-wall region are discussed 
here. With this regard, the rates of convergence of the vorticity field for three 



different window functions are presented in Fig. 13 All solutions were begun 
from the aforementioned initial conditions, and the solution were continued 
until t = 0.35. Similar time step sizes were used for all runs, corresponding 
to CFL = 0.5 for the finest grid (that is, 2048^-point grid). The first-order 
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Fig. 14. The temporal rates of convergence for the no-wah flow (left), and the flow 
under the influence of the walls (right). The rate of convergence for wall flow is 
depended on the timestep sizes. However, the maximum rate of convergence of the 
time integration method (i.e., O^df^)) was achieved. 



boundary condition setting J\fp = 1 was used, and the normalized norms 
of errors 
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were calculated in whole of the flow domain Q. As one can see, the rates of 
convergence are approached to 0{N~^) by increasing the sharpness of the 
window functions. 

In calculation of the temporal rates of convergence, and in order to reduce the 
Gibbs oscillations, a rather smooth window function f?i2(x) was used on a 512^- 
point grid, and different timestep sizes. At first, we consider the temporal rate 
of convergence of the no-wall flow. With this regard, we examine the predicted 
flow at t = 0.1, for which the effects of the solid walls are approximately 



negligible. The results are presented in the left part of Fig. [14], which shows 
the maximum rate of convergence of the fourth-order Runge-Kutta method. 
Achieving the maximum rate of convergence of the time integration method, 
primarily shows complete decomposition of the kinematics and dynamics of 
the flow in the vorticity-velocity formulation of the NSE (as it was mentioned 
in §§. 

Then we repeated the calculations for t = 0.4, where the flow is under the 
influence of the solid wall boundary conditions. The results are presented in 
the right side of Fig. [T4j With regard to this figure, the following points are 
noticeable: 

(i) Our calculations have not shown a constant rate of convergence, rather. 
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the rate of decaying increased by decreasing the timestep sizes. Similar be- 
havior were observed in our other numerical experiments on different test 
cases and a fairly wide range of timestep sizes. It seems that this behav- 
ior of the errors is a consequence of the method of implementation of the 
boundary conditions, and our insistence on the explicit formulation. More or 
less similar non-constant temporal rates of convergence have been obtained 
previously in some other embedded boundary methods [2T1[T|2] . Although 
absence of a constant rate of convergence is a drawback for a method in 
general, achievement of the highest rates of convergence of the time inte- 
gration method (that is, the fourth-order), can be seen as a merit of the 
method. 

(ii) Although for the sufficiently smooth solution fields the fourth order rate 
of convergence was obtained, our other numerical experiments have shown 
that for non-smooth solutions (obtained by e.g. sharper windowing), the 
Gibbs oscillations limits again the temporal rates of convergence to (9 (At)"* 
were 1.1 < m < 1.3. Filtering of the solution yields higher rates which are 
depended mainly on the order of the filter. 



3.1.4 Performance of the method 

The computational costs of the present method in comparison to the classical 
pseudo-spectral method is given in Tab. [2] In this table, the data are presented 
per boundary condition setting (each Runge-Kutta sub-step). Note that the 
de-aliasing process (which means zero-padding), is just needed in obtaining 
Ni and N2 in box (PS) in Fig. |2} Therefore, the FFTs in the boundary 
condition settings are performed on the original non-padded iV-point grid, not 
on the padded grid (with the benifit of increasing the efficiency of the method). 
To highlight the issue, the padded and non-padded FFTs are separated in 
Tab. [2] As one can see, implementation of the immersed boundaries needs 
just four non-padded FFTs, for each sub-step. The CPU times and required 
memories are obtained for a 3GHz Intel Pentium 4 system. Both methods used 
the same FFT subroutine and optimization levels, and the CPU times are in 

Table 2 

Computational costs for the classical pseudo-spectral (CPS), and the proposed 
Fourier immersed boundary method (FIBM); for each sub-step of fourth-order 
Runge-Kutta method. The FFTs are divided into the FFTs on the padded domain 
(P), and the non-padded domain (N-P). 



Method 


N 


#FFT (P) 


#FFT (N-P) 


CPU time 


Memory 


FIBM 


512^ 


16 


4 


0.947 Sec. 


105 MB 


CPS 


5122 


16 





0.905 Sec. 


105 MB 
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Fig. 15. Geometric parameters of the oscillating circular cylinder problem 



average. Apparently, the method is really efficient with less than 5% additional 
CPU time and without any sensible increasing in the required memory, in 
comparison to the Fourier pseudo-spectral method. Such a performance, in 
addition to robustness of the windowing procedure, makes the method ideal 
for the moving boundary, multi-object. . .problems, which some examples can 
be seen in the sequel. 



3.2 Oscillating circular cylinder in fluid at rest 



Undoubtedly, one of the most attractive features of the immersed boundary 
methods is the possibility of simple and efficient treatment of the complex and 
moving boundary problems. In the present section, to show the capability of 
the method in simulation of the moving boundary flow, the problem of in-line 
oscillating cylinder in a quiescent fluid is examined. 



3.2.1 The problem setup 

The problem has been analyzed several times experimentally [22], as well as 
numerically [3CT]l46|l33j . in a fairly wide range of involved physical parameters; 
which are the Reynolds number Re = UDu'^, and the Keulegan-Carpenter 
number KC = U{fD)^^. In these definitions, U is the maximum velocity, / 
is the oscillations frequency, and D is the cylinder diameter. In the present 
study, to facilitate assessment of the results, the well-documented combination 
(Re, KC) = (100, 5) is considered, which the previous experiments have shown 
that two-dimensional calculations are rather justifiable [22] . 
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Table 3 

The physical as well as numerical parameters, used in the solution of the oscillating 
cylinder problem. 



D 


/ 


A u 6/Ax 


N 


At 


CFL 


0.35 


1.0 


0.27852 6.1249747 x 10"^ 8 


512^ 


10-4 


0.2 



A simple harmonic oscillation in the xi direction is exerted to the cylinder 
center as 

xic(t) = Asin(27r/t), (21) 



where A is the amplitude of the oscillations. It can easily be shown that for 
such harmonic oscillations KC = 2ttAD~^. Fig. 15 illustrates the problem 
setup, and Table ([s]) summarizes the main physical and numerical parameters 
of the solution. 

The problem was solved on a fixed 512^-point equi-spaced Cartesian grid, 
distributed on a (27r)^ regular domain, using a fixed timestep size At = 2 x 
IQ-^^Sec, equivalent to CFL ^ 0.2. The zero velocity margin 6 = 8Ax was 
chosen, which resulted in 504^ active grid points. On such a grid, there were 
about 28 grid points across the cylinder diameter, and the number of numerical 
boundary points was about 30 in average. 

To keep the solution as smooth as possible, the solution was started from zero 
velocity at 6^ = — | phase angle, and therefore, a 0.25 Sec. time shifting is 
used when comparing our results with the other data. Although the no-slip 
conditions are not satisfied accurately in the initial time steps (as discussed in 



2.2, and also see [iTll^ ). it developed appropriately in the next times, and 



therefore the results for J\fj. = 1 are presented here. 



3.2.2 The captured dynamics 

The problem was solved for three different orders of boundary conditions {i.e., 
Afp = 0, 1,2), and the effects of Afp will be discussed later. Before that, the 
captured dynamics for A/^ = 2 is presented here. 



The obtained vorticity dynamics is presented in Fig. [T6j As one can see, it 
is in good agreement with the reported experimental and numerical studies 
of this (Re, KC) combination (c.f. [2^46|33j ). As the cylinder moves to one 
side, upper and lower boundary layers develop with symmetric starting and 
separation points; while a wake develops in the downstream of the cylinder, 
containing two symmetric counter-rotating vortices. Then in the opposite side 
cycle similar structures are observable, and the last wake vortices (now located 
in front of the cylinder), split by the cylinder. More quantitative comparisons 
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Fig. 16. Vorticity iso-lines of four phase angles (a): 0°, (b): 96°, (c): 192°, (d): 
The Helmholtz filter (with Ca = 0.46) is used as a post-processor in order to remove 
the Gibbs oscillations (see 



2.4). The dashed lines denote negative values. 



can be made by comparing the instantaneous velocity profiles, as illustrated 



in Fig. 17 The velocities are compared with the experimental data of Dautch 
et al. [22], and similar to some other works P^46f33j . four X\jD sections are 
illustrated for each of the three phase angles. As one can see, there is a general 
agreement between the present results and the experimental data. 
It should be noted that because of high levels of oscillations (especially for the 
vorticity fields), the aforementioned Helmholtz filter (see § 2.4) with Cq, = 
0.46 was used (as a post-processor), in obtaining the vorticities of Fig. 16 



The Ca was chosen by observation, such that the active contour levels can 
be compared with the benchmark data. Obtaining fairly acceptable results 
from an oscillating solution by use of a simple filter (just as a post-processor) 
is interesting, and it may be interpreted as a confirmation of the previous 
observations (e.g'., P7]l28f31j ) that the higher degrees of pointwise accuracies 
can be recovered from a contaminated spectral solution (at least in some 
circumstances). 



3.2.3 On the effects of order of the boundary condition setting 

Presence of various numerical and experimental studies of this problem offers 
a unique opportunity for direct evaluation of the effects of order of imple- 
mentation of the boundary conditions {i.e., A/^), on the overall accuracy of 
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Fig. 17. Comparison of the instantanious velocity profiles for three phase angles 
(a): 180° , (b): 210° , and (c): 330°. Lines show the results of the present method 
(obtained from the second-order boundary conditions Mp = 2), while the symbols 
are the experimental data |22j . 



the solution. In the last test case (z.e., the dipole-wall collision), the £^ con- 
vergence of the vorticity fields were considered, which were highly under the 
influences of the Gibbs oscillations (see ^3.1.3, and Figs. 13 and 14). Instead, 
in the present section, it is aimed to evaluate the accuracy of the solution by 
direct comparison of the velocities with a second-order solution, which has 
been verified several times (that is, the Dutsch et al. solution (22]). 
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Fig. 18. Longitudinal velocity profiles at xi/ D = of the 0° phase angle. The results 
of different Np are compared with the second-order finite volume solution of Dutsch 
et al. [22]. 



In Fig. 18 the longitudinal velocity for different orders of implementation of 
the boundary conditions (that is, Aj,=0,l,2), are compared with the second- 
order finite volume solution of Dutsch et al. [22j in the x/D = section of 
the 0° phase angle. This position is chosen intentionally because the effects 
of accuracy of the longitudinal velocity boundary condition can be evaluated 
directly. The following points can be made with regard to this figure: 

i) As one can see, the first-order and second-order implementation of the 
boundary conditions J\fp = 1, 2, yielded velocity profiles that are very close 
to the second-order finite volume solution. Therefore, we can say that for 
these orders of boundary conditions, the solution has been practically 
equivalent to a second-order solution. Comparisons of the velocities in 
the other phase angles and other x/D sections (are not presented here 
for the sake of brevity), do confirm the above statement. 

ii) Even for the least order of boundary conditions (that is Afp = 0), except 
for the regions near to the cylinder surface, the velocities are again really 
close to the second-order finite volume solutions. This issue is practically 
important, since it justifies use of easy-to-implement case J\fp = (which 
bypasses some tedious geometric operations), in rough solutions for even 
very complex geometries. 

iii) Our numerical experiments have shown that the minimum amounts of 
the Gibbs oscillations occurs for Ap = 0, and the maximum occurs for 
Afp = 2. Note that the velocities are presented here without any kind of 



filtering. In fact, just for presentation of the vorticity fields (in Fig. 16) 
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Fig. 19. Geometry of the lid-driven cavity problem. The side chanel C is considered 
in order to balance the mass flow rate and acheiving zero mean velocity and vorticity. 



the Helmholtz filter is used, because of high levels of the Gibbs oscilla- 
tions. 



3.3 Lid- driven cavity flow 



The method is not restricted to the no-shp and no-penetration immersed 
boundary conditions. In fact, any arbitrary immersed velocity boundary con- 
dition (which can be presented by a window or mask function), can be imple- 
mented via the proposed algorithm. To show this capabihty of the method, the 
classical two-dimensional lid-driven cavity flow is investigated in two different 
regimes. For the steady solutions, a fairly low Reynolds number flow Re=100 
is chosen, while the unsteady solutions is examined by analysis of a higher 
Reynolds number flow, that is Re=1000. Table |4] presents the physical and 
numerical characteristics of the solutions. 



Our suggested configuration is illustrated in Fig. [T9j As one can see, the solid 
walls of cavity are implemented by a U-shaped solid immersed body, and a 
side channel (denoted by C), is considered in order to balance the overall mass 
flow rate, and achieving zero mean velocity and vorticity fields. Both no-slip 
conditions and given velocity U are set at once in the beginning of each time 



step, as explained in § 2^, Since the cavity walls were coinciding with the 
Cartesian grid, the numerical and physical boundary points are coinciding, 
and therefore, extrapolation process is bypassed (similar to the dipole-solid 
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Fig. 20. Left: The vorticity iso-lines of the steady solution (Re=100) in comparison 
to the results of Auteri et al. [6| . The dashed lines show the results of the present 
work. Right: The velocity profiles in comparison to the results of Ghia et al. |26j 
and a classical second-order finite difference vorticity-velocity solution. 

wall problem). 

At first we examine the steady solution. Several numerical and experimental 
studies have demonstrated presence of steady solution for Re = UL/u = 100 
(see e.g. PSIISIISIIHB] ). In the present work, the problem was solved on a 2562- 
point grid, which because of presence of the margin B, the side channel C, and 
the solid U-shaped body, there were ISS^ active grid points inside the cavity. 
Since the computational cost of the method scales by (A^log A^), one can said 
that about 48% extra cost was paid for a ISS^-point simulation. The solution 
begun from zero velocity field, and a constant time step At = 2 x 10~^ Sec. 
(equivalent to CFL ^ 0.5) was used; and the solution was continued until the 
steady solution was established after about 8.6 Sec. 

On the other hand, the numerical experiments showed that Mr > 1 was 
needed for adequate satisfaction of the no-slip conditions. It can be a con- 
sequence of dominance of the solid walls on the solution in this particular 
flow. In fact, when implementing the boundary conditions, a huge amount of 
discontinuities are imposed to the u^*^ (in comparison to e.g., the oscillating 
cylinder problem). Therefore, removing these discontinuities, and enforcing the 
solenoidality, changes the immersed boundary conditions of Ug^j such that it 

Table 4 

Physical and numerical parameters of the cavity flow. 





Re 


Grid resolution 


Active grid resolution 


At 


CFL 


steady 


100 


2562 


1852 


2 X IQ-^ 


« 0.5 


unsteady 


1000 


5122 


3752 


10-4 


« 0.5 
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Fig. 21. The vorticity iso-lines of the unsteady solution for Re=1000, at the time 
instant t = 6.25 Sec, compared with the data of Auteri et al. [H]. The dashed hues 
show the results of the present method. 



is needed to set them more than once (see § 2.2). However, our experiments 
did not show changes in the results for A/'r > 3, and therefore, the present 
results are obtained with J\fr = 3. 

The results are shown in Fig. 20 In the left panel, the vorticity iso-lines are 
compared with the results of singularity-removed Chebyshev solution of Au- 
teri et al. [6\. As one can see, there is a good overall agreement between the 
results, although some discrepancies (especially in the middle of the cavity) 
are noticeable. However it should be emphasized that there are substantial 
differences between our solution and the Auteri et al. method [B], in imple- 
mentation of the boundary conditions and levels of solenoidality. In the right 
side of Fig. 20, the velocity profiles in the middle of cavity are presented and 
compared with the data of Ghia et al. [26] , and a classical second-order finite 
difference solution of vorticity-velocity formulation of the NSE. In the finite 
difference solution, the no-slip conditions are implemented via the explicit 
Thom's rule [ID], which satisfies just approximately the solenoidality [39]. As 
it can be seen, while the Ghia et al. data (which obtained from solution of 
the primitive variables formulation of the NSE, and pure Neumann boundary 
conditions for the Poisson pressure equation), are coinciding with the finite 
difference results; they show some discrepancies with our data, especially far 
from the solid walls. 
In addition to the steady solution, it is aimed to verify time accuracy of 
the method. To this end, a higher Reynolds number (that is, Re=1000) in 
an unsteady regime is examined. The unsteady cavity flow has studied sev- 
eral times experimentally and numerically (see e.g. [HIEES])- The problem was 
solved on a 512^-point equi-spaced Cartesian grid which resulted in 375^ active 
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grid points inside of the cavity. Therefore, the extra cost of the calculations 
has been about 46%. To keep the CFL number about 0.5, a constant time 
step At = 10~^ was chosen. The solution was begun from a zero velocity 
field at t = 0, and continued until t = 6.25, which the results are comparable 
by the other results [6|l5] . Like the last test case, Afr = 3 was used. In Fig. 



21 the vorticity iso-lines in comparison to the results of singularity-removed 



Chebyshev solution [6J are shown. Good agreements between the results are 
observable, despite the fact that a different time integration method (that is, 
the second-order Adams-Bashforh method) is used in obtaining the results of 
P and [5]. 



4 Conclusions 



An immersed boundary Fourier pseudo-spectral method was proposed for 
the vorticity- velocity formulation of the two-dimensional incompressible NSE. 
The zero-mean Fourier pseudo-spectral solution are used, and therefore, the 
method is applicable to the confined flows, which are modeled by considering a 
zero velocity margin in the vicinity of the regular boundaries. Without explicit 
addition of external forcing functions, arbitrary Dirichlet velocity boundary 
conditions are implemented by direct modification of the diffusion and convec- 
tion terms of the vorticity transport equation; and in this way, it was shown 
that the obtained velocities are solenoidal. The immersed boundary conditions 
are approximated on some regular grid points (called the numerical boundary 
points), by use of different orders (up to second-order) polynomial extrapola- 
tions, in the normal directions to the immersed surfaces. 
Although because of presence of the Gibbs phenomenon the spatial rates 
of convergence are not so high (even less than two, especially in the regions 
near the solid walls), good agreements between our results with some other 
second-order solutions were observed in some test cases. On the other hand, 
the temporal rate of convergence was found to be a function of timestep sizes, 
and for sufficiently small timestep sizes, the highest order of time discretiza- 
tion (that is, the fourth order) was obtained. 

Like many other immersed boundary methods, it was observed that imple- 
mentation of the immersed velocity boundary conditions imposes some dis- 
continuities; and conversely, imposing the continuity changes the immersed 
boundary conditions. However, it was observed that by repeating the proce- 
dure, and by development of the solution, some solenoidal velocities, which 
satisfy the immersed boundary conditions (approximately) are obtainable. 
The method was implemented to a conventional Fourier pseudo-spectral solver 
of the vorticity-velocity form of the NSE, without any changes in the inter- 
nal structure of the solver, and the performance of the method was compared 
with the original solver. In this context, it was observed that fixed and moving 
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immersed boundaries can be treated by a computational cost, which scales by 

{N logiV). 

The best accuracies (and also the most Gibbs oscillations), were obtained 
for the second-order boundary conditions. However, even for the easy-to- 
implement zero-order boundary conditions, fairly accurate solutions were ob- 
tained, at least far from the solid walls. Particularly this feature of the method 
is interesting, because this simple algorithm can be used in a rather wide 
range of applications, from rough simulations to the fairly accurate ones, just 
by changing the order of implementation of the boundary conditions. In our 
opinion, it is particularly a consequence of using a pseudo-spectral solver as 
the core of the method. 

Just a simple constant-width Hclmholtz filter was used in the present work, 
as a postprocessor, and in order to remove the Gibbs oscillations in visualiza- 
tion of the results. However, like other spectral methods, one can use some 
particular numerical filters to improve the rates of convergence. This is one of 
several lines which we can assume for the extension of the present method. 
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Appendix: Some fundamental considerations about the proposed 
algorithm 

Since the zero-mean Fourier series are used in the present method, the method 
is applicable to the flows with zero-mean velocity and vorticity, and with zero 
dynamics of the mean vorticity (see Eqns. Q or (|8|). Now from the Stokes 
theorem, the mean vorticity reads 



a;(0,0) = 




Among many flow conflgurations with a) = 0, the ones that ur^ = are dealt 
in the present work, and these flows are modeled by considering a zero- velocity 
margin, as it is shown in Fig. [T] 

With this in mind, the present appendix is devoted to explanation of three 
fundamental issues, particularly related to steps 3 and 4 of the proposed al- 



gorithm (§ 2.2). These discussions are presented here separately for the sake 



of clarity, and to prevent disturbing the main flow of the paper. 



.1 The mean value of the conditioned vorticity 



Clearly, Eq. (Rl yields a zero-mean conditioned vorticity (since ki = k2 = 
results in u^'^O, 0) = 0), as it is desired. However, legitimacy of application of 
this equation for arbitrary u"^*" (obtained from arbitrary immersed boundary 
conditions and different window functions), needs some considerations. 
In general, the Stokes theorem states that for a confined flow (ur^ = 0) the 
mean vorticity is vanished, regardless of distribution of velocities in D (see 



Eq. .1). However, the point is that the Stokes theorem is applicable to the 
velocity fields that are at least; the condition which can be violated in the 
general immersed boundary condition settings and windowings. 
On the other hand, note that it is desired to find the solution in the Fourier 
space, and regardless of rate of decaying of u^*" = FT{u^*-'} (which is de- 
pended on the smoothness of u^*-^), however, FT~^{u^^} has C°° smoothness 
(since it is constructed from the Fourier functions which are C°° smooth). 
Now, by applying the Stokes theorem on FT~^{FT{u^^}}, one can find 

^ I (V X n'''')dQo = ^ / FT-i{FT{u^C||^r^)^r^_ 



On the other hand, it should be noted that FTju^*-'} can be suffered from the 
Gibbs oscillations because of presence of discontinuities in u^^, and therefore, 
FT"^{FT{u^^}} is not necessarily zero on F/). Consequently, is not van- 
ishing in general. 
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However, since the zero-mean flows are dealing, these supposed to be 
discrepancies from the flow physics, and they are ignored in the solution pro- 
cedure. On the other hand, as our numerical experiments were shown, these 
spurious mean vorticities decrease as the solution is developed; and for the log- 
ical timestep sizes, they get practically-ignorable values almost immediately 
after beginning of the solution. 



.2 On the solenoidality of the velocities 



A key feature of the present method is that the physical (solenoidal) velocities 
are obtained and used, regardless of complexity of the immersed boundary 
conditions. Here it is desired to justify why Eq. ([T]) results in solenoidal veloc- 
ities for an arbitrary cu^*^. 



Proposition 1 For any conditioned vorticity u^*^ G C?{D) the velocity vector 
uio^i = -^t!^^"" (-3) 



k 



is solenoidal, in which Co^'^ = FTju;^^}. 

Note that Eq. Q is re-written again (that is Eq. [3]) for clarity. 

PROOF. There are many ways to proof the above proposition, and we give 

one of them which is more consistent with our purposes. 

At first, since cu^^ e C^{D), the Fourier coefficients a)^^(k) = FT{a;^^} exist 

for any k, although their rate of decaying can be not so high. Now, it is a 

well-known fact that the continuity equation in the Fourier space reduces 

to 

k ■ usoi = 0. (.4) 
In the other words, the solenoidal velocities are perpendicular to the wavenum- 



ber vector and vice versa. Now referring to Eq. (.3), the proposition is proven 
since the velocities Ug^j are obtained from k-*-, and therefore they are perpen- 
dicular to k. 
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In addition to the above proof, also our extensive numerical experiments on 
different test cases have shown that the divergence of u^^i remains in the order 
of machine accuracy everywhere in D. 

The above simple proof in the Fourier space should not hide the detailed 
mechanism that discards the divergence. In fact, by taking the divergence of 
both sides of Eq. ^ we have 



which means the obtained velocities from Poisson equation Eq. ([2]) are such 
that the divergence of them is a harmonic function (and therefore gets its max- 
imum values on the boundary r^), regardless of the boundary conditions. As 
an immediate consequence, the divergence of the velocities can be controlled 
by controlling the divergence on (near) the boundary Td. This statement is 
not new, and it has been mentioned and discussed previously in many refer- 
ences (see e.g. |20] or [39]). 

Now definition of a zero velocity margin B in the vicinity of To can be seen 
as a way for controlling the divergence. This strategy is not restricted to the 
spectral method, and can be used in other numerical methods, e.g. the finite 
difference method. However, in the Fourier pseudo-spectral solutions, since 
there is not any boundaries, the divergence vanishes perfectly everywhere on 
D, as it proved. In these cases the margin B helps achieving the higher rates 
of decaying of the Fourier coefficients. 

In fact, we found that considering a closed flat zero velocity margin in the 
vicinity of the regular boundary (when it is possible), is a good idea for over- 
coming many difficulties related to finding appropriate vorticity boundary 
conditions and obtaining the solenoidal velocities. We will follow this idea in 
our next works on the fluid-solid interaction problems. 



.3 About the dynamics of the mean vorticity 

There is still another question that should be answered with regard to the 
proposed algorithm. In fact, by time integration of Eq. ([s]), we assumed zero 
dynamics for the mean vorticity during the time interval At. Legitimacy of 
this assumption is in order. 

Note that in modeling of the flow conflguration of Fig. [T| the modifled L, 
Ni, and N2 are introduced into the vorticity transport equation, and then 
this modifled equation is transformed into the Fourier space (yielded Eq. (|8|). 
Therefore, in order to analysis of the dynamics of the mean vorticity of our 
numerical model, the mean mode of the modifled vorticity transport equation 
should be analyzed. 

Theorem 1 Given a conditioned vorticity o;^'" = V x u^*-^, time integration 





36 



results in zero dynamics for the mean vorticity, in which u|^[ = FT ""^{usoi}' 
and ufoi is obtained from Eq. 

PROOF. To show dtOJ = 0, one way is to show that the modified L, Ni, and 
N2 have zero means. Now, at first, note that the diffusion term L (that is, 
the Laplacian of vorticity) is essentially silent about the mean values, since 
it consists of the (second-order) derivatives in both directions. On the other 
hand, the mean values of the other two convection terms (that is, Ni, and 
N2) can be found directly. For example 

(Ni)io^. = ^V-^'V-^\N,)lg = ^V-^V-^'[V^,V^,{ul - uDlg], (.7) 



where stands for derivative with respect to x^; while "D^.^ = Jq dxi shows 
the inverse of derivative operator (that is, integration on the regular do- 
main D). Now, implementation of the integral operators and changing the 
order yields 



(Ni)io^. = 7^(^^-^2^.J[2^;,^^^x.(«^ -«?)io'^i] 

= 7V[(«2-«.l!)^]l^ (.8) 

On the other hand, since u|^j is double periodic on Z), one can write 

{ul - Mi)ioi(a;i = 0, X2) = {ul - Mi)|^i(a;i = h, xa), 

(«2 - w?)ioi(a;i>a;2 = 0) = {ul - ul)lg{x,,X2 = h). (.9) 



Substitution of these relations into Eq. (.8), yields (Ni)|qi = 0. 

Exactly the same procedure can be followed for the (N2)soi "which results in 

(N2)ioS = 0. (.10) 
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Remark 2 Note that changing the derivatives order from Eq. (.7) to Eq. (.8) 
needs uf^j G C^{D) at least. On the other hand, as it was argued before, u|^^ 
FT"^{FT{ — V X w^'-'}} G C°°(Z)) regardless of the rate of convergence of 

-BC 
%ol- 

Remark 3 Although the theorem and proof are offered for this particular 
form of the vorticity transport equation, that is Eq. (|3]); it is worth mention- 
ing that they are valid and applicable to the conventional form of vorticity 
transport equation ([T]), but it is just needed applying the Gauss theorem once, 
on each convection term. 
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